Question 1

Identifying the Area of Interest

library(raster)
library(tidyverse)
library(getlandsat)
library(sf)
library(mapview)
library(osmdata)

bb = read_csv("../data/uscities.csv") %>%
  filter(city == "Palo") %>%
  st_as_sf(coords = c("lng", "lat"), crs = 4326) %>%
  st_transform(5070) %>%
  st_buffer(5000) %>%
  st_bbox() %>%
  st_as_sfc() %>%
  st_as_sf()

bbwgs = bb %>%
  st_transform(4326)

Question 2

Identifying, Downloading, and Caching Images

meta = read_csv("../data/palo-flood.csv")

files = lsat_scene_files(meta$download_url) %>% 
  filter(grepl(paste0("B", 1:6, ".TIF$", collapse = "|"), file)) %>% 
  arrange(file) %>% 
  pull(file)

st = sapply(files, lsat_image)

b = stack(st) %>% 
  setNames(c(paste0("band", 1:6)))

cropper = bbwgs %>% 
  st_transform(crs(b))

r = crop(b, cropper)
Dimensions, Coordinate Reference System, and Resolution of stacked images
dim(b)
## [1] 7811 7681    6
crs(b)
## CRS arguments:
##  +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs
res(b)
## [1] 30 30
dim(r)
## [1] 340 346   6
crs(r)
## CRS arguments:
##  +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs
res(r)
## [1] 30 30

Question 3

Making RGB plots of these images.

par(mfrow = c(2, 2))
plotRGB(r, r = 4, g = 3, b = 2)
plotRGB(r, r = 5, g = 4, b = 3)
plotRGB(r, r = 6, g = 5, b = 4)
plotRGB(r, r = 5, g = 7, b = 1)

par(mfrow = c(2, 2))
plotRGB(r, r = 4, g = 3, b = 2, stretch = "lin")
plotRGB(r, r = 5, g = 4, b = 3, stretch = "lin")
plotRGB(r, r = 6, g = 5, b = 4, stretch = "lin")
plotRGB(r, r = 5, g = 7, b = 1, stretch = "lin")

par(mfrow = c(2, 2))
plotRGB(r, r = 4, g = 3, b = 2, stretch = "hist")
plotRGB(r, r = 5, g = 4, b = 3, stretch = "hist")
plotRGB(r, r = 6, g = 5, b = 4, stretch = "hist")
plotRGB(r, r = 5, g = 7, b = 1, stretch = "hist")

Applying a color stretch helps make differences in landscapes more visibly clear to assist in identifying patterns or features.

Question 4

Assessments of Surface Water Features Using A Threshold Value.

Raster Algebra for 5 Formulas

ndvi = (r$band5 - r$band4) / (r$band5 + r$band4)
ndwi = (r$band3 - r$band5) / (r$band3 + r$band5)
mndwi = (r$band3 - r$band6) / (r$band3 + r$band6)
wri = (r$band3 + r$band4) / (r$band5 + r$band6)
swi = 1 / sqrt(r$band2 - r$band6)

sw_rasters = stack(ndvi, ndwi, mndwi, wri, swi) %>% 
  setNames(c("ndvi", "ndwi", "mndwi", "wri", "swi"))

palette = colorRampPalette(c("blue", "white", "red"))
plot(sw_rasters, col = palette(256))

The images are similar in that they highlight the river region as well as other surface water areas. They also all use the the same blue, white, and red color scheme that we created in our palette. They differ in that the surface water areas are not all highlighted in the same color, and some images direct more attention towards the landscape rather than water regions. The final plot, the simple water index, highlights only the surface water captured in the image.

Raster Thresholding

thresh_ndvi = function(x){ifelse(x <= 0, 1, NA)}
thresh_ndwi = function(x){ifelse(x >= 0, 1, NA)}
thresh_wri = function(x){ifelse(x >= 1, 1, NA)}
thresh_swi = function(x){ifelse(x <= 5, 1, NA)}

flood_ndvi = calc(ndvi, thresh_ndvi)
flood_ndwi = calc(ndwi, thresh_ndwi)
flood_mndwi = calc(mndwi, thresh_ndwi)
flood_wri = calc(wri, thresh_wri)
flood_swi = calc(swi, thresh_swi)

flood = stack(flood_ndvi, flood_ndwi, flood_mndwi, flood_wri, flood_swi) %>% 
  setNames(c("ndvi", "ndwi", "mndwi", "wri", "swi"))

plot(flood, col = "blue")

flood <- na.omit(flood)

Question 5

Supervised and Unsupervised Classification

set.seed(09072020)
dim(r)
## [1] 340 346   6
values = values(r)
idx = which(!is.na(values))
v = na.omit(values)
vs = scale(v)

library(stats)

k12 = kmeans(vs, centers = 12)
kmeans_raster = flood$ndvi
values(kmeans_raster) = k12$cluster
plot(kmeans_raster)

k11 = kmeans(vs, centers = 11)
kmeans_raster = flood$ndvi
values(kmeans_raster) = k11$cluster
plot(kmeans_raster)

k10 = kmeans(vs, centers = 10)
kmeans_raster = flood$ndvi
values(kmeans_raster) = k10$cluster
plot(kmeans_raster)

k13 = kmeans(vs, centers = 13)
kmeans_raster = flood$ndvi
values(kmeans_raster) = k13$cluster
plot(kmeans_raster)

k9 = kmeans(vs, centers = 9)
kmeans_raster = flood$ndvi
values(kmeans_raster) = k9$cluster
plot(kmeans_raster)

flood_values = getValues(flood_ndvi)
kmeans_values = getValues(kmeans_raster)
table = table(flood_values, kmeans_values)
which.max(table)
## [1] 2
mask = function(x){ifelse(x == 5, 1, NA)}
flood_mask = calc(kmeans_raster, mask)

flood = addLayer(flood, flood_mask)

Question 6

s1 = cellStats(flood$ndvi, stat = 'sum') * 30^2
s2 = cellStats(flood$ndwi, stat = 'sum') * 30^2
s3 = cellStats(flood$mndwi, stat = 'sum') * 30^2
s4 = cellStats(flood$swi, stat = 'sum') * 30^2
s5 = cellStats(flood$wri, stat = 'sum') * 30^2
s6 = cellStats(flood_mask, stat = 'sum') * 30^2
comparison = cbind(method = c("ndvi", "ndwi", "mndwi", "wri", "swi", "kmeans_mask"), cell_area = c(s1, s2, s3, s4, s5, s6))


knitr::kable(comparison, 
             caption = "Total Area of Flooded Cells",
             col.names = c("Band", "Total Flooded Cells")) %>% 
  kableExtra::kable_styling("striped", full_width = TRUE)
Total Area of Flooded Cells
Band Total Flooded Cells
ndvi 5999400
ndwi 6490800
mndwi 10745100
wri 13680900
swi 7622100
kmeans_mask 8017200 
total_flood = calc(flood, sum)

plot(total_flood, col = RColorBrewer::brewer.pal(9, “Blues”))

plot shows up in R but shows up as empty plot in Rmd.

total_flood <- na.omit(total_flood) flood = addLayer(flood, total_flood) mapview(total_flood)